% clear all, close all
% disp(' ');

%AM wave DSB with carrier
figure,
Fs = 75;
fc = 1; % input('Enter the carrier signal freq in hz,fc =');
fm= 0.3; %input('Enter the modulating signal freq in hz,fm =');
m= 0.1; % input('Modulation index,m= ');
n=0:1/Fs:60;
c=sin(2*pi*fc*n);%carrier signal
M=cos(2*pi*fm*n);% modulating signal
y=(1+m*M).*c;%AM signal
yd = downsample(y,5);
time = downsample(n,5);
subplot(2,1,1);plot(n,y);
ylabel('amplitude');xlabel('time index');
hold on
plot(time,yd,'r')
subplot(2,1,2); plot(n,m*M); axis([0 60 -.5 .5])
ylabel('amplitude');xlabel('time index');

% First example using a periodic cov function
% Here are the covariance contributions:
covfunc = {@covSum, {@covPeriodic, @covPeriodic}};    % a quasi-periodic component (the amplitude of the oscillations varies)

meanfunc = @meanZero;
tic
% Select priors for the covariance parameters. To evaluate the number of
% parameters in the covariance function use: feval(covfunc{:})
hyp.cov = [1 0 0 1 1 -1]; hyp.lik = -2;               % init hypers

x = time(time<=45);
x=x(:);
y = yd(time<=45);
y= y(:);
yy = yd(time>45);
yy = yy(:);
xx = time(time>45);
xx = xx(:);
% Fit the GP model
[hyp fX i] = ...     
    minimize(hyp, @gp, -500, @infExact, meanfunc, covfunc, @likGauss, x, y);

% covfunc1 = @covPeriodic; hyp1.cov = hyp.cov(1:3); hyp1.lik = hyp.lik;
% covfunc2 = @covPeriodic; hyp2.cov = hyp.cov(4:6); hyp2.lik = hyp.lik;

% Make predictions 10 years into the future
zz = n';
% hyp.cov(3) = 0;
% hyp.cov(5) = 0;
[mu s2] = gp(hyp, @infExact, [], covfunc, @likGauss, x, y, zz);
% 
% [mu1 s21] = gp(hyp1, @infExact, [], covfunc1, @likGauss, x, y, zz);
% [mu2 s22] = gp(hyp2, @infExact, [], covfunc2, @likGauss, x, y, zz);

% Plot the data and the predictions
figure; 
% subplot(311)
f = [mu+2*sqrt(s2); flipdim(mu-2*sqrt(s2),1)];
fill([zz; flipdim(zz,1)], f, [7 7 7]/8); hold on; 
plot(zz, mu,'--k'); plot(x,y,'b.'); plot(xx,yy,'r.');                            
xlabel('Year'); ylabel('CO_2 concentration, ppm');
title('GP-1 sexp+periodic+sexp+noise');
% subplot(312)
% f = [mu1+2*sqrt(s21); flipdim(mu1-2*sqrt(s21),1)];
% fill([zz; flipdim(zz,1)], f, [7 7 7]/8); hold on; 
% plot(zz, mu1,'--k');
% subplot(313)
% f = [mu2+2*sqrt(s22); flipdim(mu2-2*sqrt(s22),1)];
% fill([zz; flipdim(zz,1)], f, [7 7 7]/8); hold on; 
% plot(zz, mu2,'--k');

fprintf('\nHR = %d vs %d\n',round(60/exp(hyp.cov(2))),round(fc*60))
fprintf('\nRR = %d vs %d\n',round(60/exp(hyp.cov(5))),round(fm*60));

toc
%%

%FM generation
%fc=input('Enter the carrier signal freq in hz,fc=');
%fm=input('Enter the modulating signal freq in hz,fm =');
fc = 1; % input('Enter the carrier signal freq in hz,fc =');
fm= 0.3; %input('Enter the modulating signal freq in hz,fm =');
m=1;
t=0:0.01:60;
c=sin(2*pi*fc*t);%carrier signal
M=sin(2*pi*fm*t);% modulating signal
subplot(3,1,1);plot(t,c);
ylabel('amplitude');xlabel('time index');title('Carrier signal');
subplot(3,1,2);plot(t,M);
ylabel('amplitude');xlabel('time index');title('Modulating signal');
y=sin(2*pi*fc*t+(m.*sin(2*pi*fm*t)));
subplot(3,1,3);plot(t,y);
ylabel('amplitude');xlabel('time index');title('Frequency Modulated signal');
yd = downsample(y,5);
time = downsample(t,5);
% Here are the covariance contributions:
covfunc = @covQPeriodic;    % a quasi-periodic component (the amplitude of the oscillations varies)

% Select priors for the covariance parameters. To evaluate the number of
% parameters in the covariance function use: feval(covfunc{:})
hyp.cov = [0 .5 1 1 -1]; hyp.lik =5;               % init hypers
X = time(:);
Y = yd;
Y=Y(:);
x = time(time<=30);
x=x(:);
y = Y(time<=30);
y=y(:);
yy = Y(time>30);
yy=yy(:);
xx = time(time>30);
xx = xx(:);

% Fit the GP model
[hyp fX i] = ...     
    minimize(hyp, @gp, -50, @infExact, [], covfunc, @likGauss, X, Y);

% Make predictions 10 years into the future
zz = t';
[mu s2] = gp(hyp, @infExact, [], covfunc, @likGauss, X, Y, zz);

% Plot the data and the predictions
figure; 
f = [mu+2*sqrt(s2); flipdim(mu-2*sqrt(s2),1)] + mean(y);
fill([zz; flipdim(zz,1)], f, [7 7 7]/8); hold on; 
plot(zz, mu,'--k'); plot(x,y,'b.'); plot(xx,yy,'r.');                            
xlabel('Year'); ylabel('CO_2 concentration, ppm');
title('GP-1 sexp+periodic+sexp+noise');
